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Abstract 

We theoretically investigate spin-dependent carrier dynamics due to the electron-phonon inter- 
action after ultrafast optical excitation in ferromagnetic metals. We calculate the electron-phonon 
matrix elements including the spin-orbit interaction in the electronic wave functions and the in- 
teraction potential. Using the matrix elements in Boltzmann scattering integrals, the momentum- 
resolved carrier distributions are obtained by solving their equation of motion numerically. We find 
that the optical excitation with realistic laser intensities alone leads to a negligible magnetization 
change, and that the demagnetization due to electron-phonon interaction is mostly due to hole 
scattering. Importantly, the calculated demagnetization quenching due to this Elliot- Yafet type 
depolarization mechanism is not large enough to explain the experimentally observed result. We 
argue that the ultrafast demagnetization of ferromagnets does not occur exclusively via an Elliott- 
Yafet type process, i.e., scattering in the presence of the spin-orbit interaction, but is influenced to 
a large degree by a dynamical change of the band structure, i.e., the exchange splitting. 

PACS numbers: 75.78.Jp, 75.70.Tj, 78.47.J- 
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I. INTRODUCTION 



It was first demonstrated more than ten years ago that the magnetization of ferromag- 
nets can be "quenched" on ultrashort time scales after ultrafast optical excitation.- Apart 
from the possibilities for the ultrafast manipulation of ferromagnetism in applications, this 
observation raised the question how demagnetization dynamics in ferromagnets on a time 
scale of a few hundred femtoseconds can be understood. 

Aside from the phenomenological three-temperature model,-^ which leads to quite suc- 
cessful comparison with experiment, there are several theoretical models and experimental 
results that try to explain aspects of the underlying microscopic dynamics. For instance, 
the analysis of X-ray magnetic circular dichroism measurements suggested that the orbital 
magnetic moment does not play a prominent role in the demagnetization dynamics.-' 4 The 
authors of Ref. [3 concluded that an ultrafast spin-lattice coupling should be operative to 
explain the results. It has also been argued based on experimental results 5 that the exci- 
tation of magnons should play an important role. On the theory side, magnetic switching 
due to electronic transitions during the duration of a pump laser pulse has been analyzed 
in ferromagnets^ as well as in oxides (including phonons), 8 and the Landau-Lifshitz-Bloch 
equations have been used to describe the magnetic dynamics.- 

Perhaps the most popular microscopic explanations of the effect involve variations of 
the so-called Elliott- Yafet mechanism, in which demagnetization (or depolarization in semi- 
conductors) is due to incoherent scattering of carriers between states that are spin-mixed 
due to the spin-orbit interaction. Electron-electron scattering is a possible candidate as 
the underlying scattering mechanism,— but the focus is usually on the effects of electron- 
phonon scattering in quasi-equilibrium.— ~— In addition, super diffusive transport processes 
can contribute to the measured Kerr-effect signal because minority and majority electrons 
may simply leave the probe area at different speeds.^- 5 " 

In our opinion, it has so far been impossible rule out a single one of these mechanisms, let 
alone to pinpoint the dominant one. As a first step in this direction, we analyze a parameter- 
free microscopic model for ultrafast demagnetization and compare it with experiment. To 
keep things simple while allowing a conclusive statement, we exclusively treat the effect of 
the optical excitation and electron-phonon scattering at the level of Boltzmann scattering 
integrals while neglecting dynamical changes in the band structure, i.e., the exchange split- 



ting, in the course of the dynamics. We evaluate the model using ab-initio results for the 
simple ferromagnets nickel and iron. Comparison of the results of the present paper with 
experimental data, which are available from many different measurements, will show that a 
model without band structure changes yields a demagnetization that is too small. 

This paper is organized as follows. In Sec. II we present the dynamical equations for 
the carrier distribution functions and show how we calculate the electron-phonon and dipole 
matrix elements using a first-principles approach. In Sec. Ill we discuss the numerical 
results of this model and show that the demagnetization using realistic parameters for the 
ultrashort-pulse excitation is due to hole dynamics, but too small to agree with experiment. 
A qualitative consideration shows that this conclusion should not be altered by including 
additional scattering processes. Sec. IV contains the conclusions and the Appendix describes 
details of our numerical evaluation of the dynamical equations. 



II. THEORY 



A. Dynamical equations 

The basic idea of this paper is to integrate the dynamical equations of motions for the 
band- and momentum- resolved carrier distributions n'tit). Our model includes the inco- 
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herent scattering dynamics due to the electron-phonon interaction, as well as the optical 
excitation, so that the general form of the dynamical equation is 
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The optical excitation is given by 
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Here, is the energy of a carrier in a single-particle state tp'i with band index ji and 
momentum k. The dipole matrix element for a transition connecting two such states is 
denoted by dt u = ef iV'jj*)- The optical excitation is characterized by the dynamical 
electric field amplitude E(t), a central laser frequency Co>l, and the function g(e) that includes 
the spectral profile of the laser pulse. 

The electron-phonon contribution to the carrier dynamics at the level of Boltzmann 
scattering integrals reads 
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have contributions from absorption ("+") and emission ("— ") processes. The (angular) 
frequency of a phonon mode A is designated by u>~ and its occupation at quasi-momentum q 
by rig. The electron-phonon interaction matrix elements M± X (k', v — > k,fx) result from the 
change of the electron-lattice interaction energy due to the vibrational motion of the nuclei. 
For small displacements they are given by 
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The electron-lattice interaction potential V(r; {Ri}) depends on the electron position r and, 
in principle, on the set of the positions of the nuclei {Ri} in the crystal composed of N unit 
cells with atomic mass M. The polarization vector of the phonon mode (q, A) is denoted by 
7?^. The upper and lower signs in the exponential are associated with the phonon absorption 
M+ X (k', v — > k,fi) and emission M^ x (k', v — > k,[i) matrix elements, respectively. 

In our calculation, we assume that the phonon occupation numbers are time independent 
and remain at their equilibrium values 

This amounts to a bath assumption for the phonon system, and in this paper we fix its 
temperature at T = 300 K, i.e., the temperature of the unexcited system in most studies of 
demagnetization dynamics. 

From the dynamical electronic occupation numbers calculated according to Eq. (TTJ, we 
obtain the time-dependent magnetization by 

M(f) = ^££W2 n ?(*)' (7) 



with the single-particle spin expectation value 

wj= W*2> (8) 

in the Bloch state and the Bohr magneton [i^- In writing these relations, we have chosen 
the z-direction as direction of the ferromagnetic polarization. Orbital angular-momentum 
contributions to the magnetization are neglected as their influence on the magnetization of 
elementary ferromagnets is small. 



B. Electron-phonon matrix elements 

The numerical evaluation of Eq. ([[]) requires as input material properties, in particu- 
lar, the electronic band structure et, the spin expectation value of the single-particle states 
(S z )p the dipole transition matrix elements dt u , the phonon dispersion and, impor- 
tantly, the electron-phonon matrix elements M± (k', v — > k, /J,). We obtain these quantities 
from density- functional theory to avoid the introduction of adjustable parameters. To this 
end, we employ the augmented spherical wave (ASW) method 16 as described in the mono- 
id). The implementation of the ASW method used by us was 
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17 (see also Ref. 



developed in the Kiibler group and relies on the scalar relativistic and local spin-density 
approximations. It includes spin-orbit coupling in a second variational correction. 

The starting point for the calculation of the matrix elements is the representation of the 
wave function of a single-particle state with band index /i and momentum k in the ASW 
basis. Since we employ the atomic sphere approximation (ASA) it is usually sufficient to 
know the wave functions inside the atomic spheres where they are given by 

^if) = Y,[ C ^)i l hUr) + At(k)i l jia(r)}Y L (r) Xc r (9) 

L,cr 

where hi a (r) and ji a (r) are augmented spherical waves, Yl{t) spherical harmonics, and x<? 
Pauli spinors. Here, L = (I, m) is a multi-index that includes both the angular momentum 
and the magnetic quantum number. The relatively simple expression given here is only 
valid for materials with basis consisting of a single atom, such as the simple ferromagnets 
investigated in the present paper. The spherical wave functions, together with the coefficients 
A^ La {k) and C£ a (k), are calculated self-consistently during the iterative solution of the Kohn- 
Sham equations. 



For the evaluation of Eq. (jSJ), we employ the so-called rigid ion approximation. That is, 
we assume that we can write the lattice-configuration dependence of the electron-phonon 
interaction potential 

V{f-{R l }) = Y J <r-R l ) (10) 
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as a superposition of the on-site potentials v(r), see below, Eq. ( TT2l) . We also assume that 
the potential v vanishes outside the atomic sphere. The rigid ion approximation is known 
to give a quite realistic description of the electron-phonon coupling in transition metals.— 
Equation ([5]) can then be simplified to yield 

Mf A (£>^ k,n) = - J^^xh^-k,G J vc d 3 r Vf (r) ■ V^r)]^,^, (11) 

where UC denotes an integral over the unit cell, which due to the ASA is assumed to 
be spherical. For the on-site potential experienced by the electrons, we include the spin- 
averaged radial Kohn-Sham potential V e s(r) as well as the spin-orbit interaction 

v(r) = V cS (r) + — — 2-— — - a-(rxp). (12) 
{2mc) t ar 

The additional spin-orbit term is often neglected for the electron-phonon interaction, even 
though it has been shown to be of importance for spin-relaxation in materials with time- 
inversion symmetry— Not much is known about the influence of this term in ferromagnets 
where the time-inversion symmetry is broken. We therefore calculate the matrix element 
with and without the spin-orbit contribution and show that our final results on ultrafast 
demagnetization are not qualitatively influenced by the inclusion of the second term. 

We first give the result for the calculation without the spin-orbit contribution in Eq. ( fl"2l) . 
In this case we can directly evaluate the integral over the unit cell in Eq. (jTTT) : 
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The radial matrix elements 



(f\^\9)i^ = (-iyi W I r 2 f la {r)^^~g Va ,{r)dr (14) 



can be calculated by integrating the gradient of the Kohn-Sham potential and 

G w = [ dnY*(r)rY Lf (r) (15) 



can be evaluated in terms of Gaunt coefficients.— 

The calculation of the electron-phonon interaction matrix element flTT]) including the 
spin-orbit contribution could, in principle, be achieved by evaluating the integral 
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and adding it to Eq. ( fl3l) . However, straightforward numerical evaluation of Eq. ( f!6l) runs 
into difficulties because of the strong divergence of the integrand for r — > 0. We circumvent 
this problem by calculating the complete matrix element (fTTj) by rewriting the gradient of 
the potential including the spin-orbit interaction [Eq. (fl2l) ] as a commutator 

Vv(7)= l -\p,H eS ] (17) 
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with the Hamiltonian H e g = 2^ + v(r). For the evaluation of matrix elements of H e s, we will 
assume that it produces the energy eigenvalues when acting on the corresponding Kohn- 
Sham eigenvector, even though the eigenenergies and eigenvectors are computed using scalar 
relativistic corrections to H e g. If we now assume completeness of the ASW basis, Eq. ([Tj 
may be used for a reformulation in terms of the momentum matrix elements 



d 3 rM* (r)Vu(r)V£(^ = 4 / d 3 r ^*(f) [pH cS - H eS p\ ^,{f) 

..<■ ^ ^ h Ji:c ^ (18) 



where the momentum matrix elements (V^qPlV^uc are ca l cu l a tecl in the ASA using a 
consistent method developed by Oppeneer et al.— Since Vv(f) is a hermitian operator on 
the unit cell, we can derive the following expression that is more symmetric with regard to 



initial and final states: 
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Although our assumption of completeness may yield the matrix element only with a certain 
error because the ASW method uses a rather small number of basis functions, the qualitative 
conclusions discussed in the next section are not affected by this. 

In the numerical calculations, we typically used a fc-point grid of about 2000 points in 
the irreducible wedge of the band structure, and the dynamical equations were solved on the 
same grid (see the Appendix n for details on the numerical method). Experimental values for 
the lattice constants^ were used. The phonon dispersion was calculated with Quantum 
Espresso^ in the same way as by Dal Corso et al.— The latter paper shows that phonon 
dispersions obtained in this approach are in good agreement with experimental data. 



C. Dipole matrix-elements 

The dipole matrix elements are calculated by reformulating them in terms of the momen- 
tum matrix elements 



S T = = ^ (9 x W eff M) |4 (20) 
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where the momentum matrix elements (ipt \p\^) are again calculated according to Oppeneer 
et al.— The contribution of the spin-orbit interaction is usually neglected. We include it here 
as it may directly contribute to spin flips, even though our numerical results in the end will 
show that the difference is insignificant. It can be calculated from the wavefunctions as 
follows: 



h 



8 



V£ x Weff(r)) 



V4 



X ^LL' 



o\ct' L,L' 



d 3 r^f(r)(axf) ^%(r) 



(21) 



cf. Eq. (I14p for the definition of the matrix elements involving V e g. 
III. RESULTS 

In this section, we present numerical results obtained from the solution of the dynamical 
equation ([I]) and some qualitative considerations on the role of scattering processes in ultra- 
fast demagnetization dynamics. For the calculations we use matrix elements computed as 
described in the previous sections. Further details of our numerical implementation of the 
dynamics are included in the Appendix. 

A. Optical excitation 



We first examine the excitation process by the ultrashort optical pulse without including 
the electron-phonon interaction. We model a homogeneous excitation of a ferromagnetic 
metal by a laser pulse with gaussian temporal shape, full width at half maximum of 50 fs, and 
a spectral width of 100 meV at a central photon energy of ftw = 1.55 eV. These parameters, as 
well as the pulse intensity of 4mJ/cm 2 , are chosen to match typical experimental excitation 
conditions. To determine the electric field amplitude present in the material we note that 
the chosen intensity corresponds to an electric field amplitude of Eq = 7.5 x 10 8 V/m in 
vacuum. Reflection at the surface as well as the optical density of the material lead to a 
reduction of the field amplitude in the material 
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En = E(] 



(22) 



where n and k denote the real and the imaginary part of the refractive index, respectively. 
Taking n = 2.22, k = 4.90 for nickel and n = 2.92, k = 3.36 for iron, this leaves us 
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Figure 1. (Color online) Magnetization dynamics due to optical excitation alone. The magnetiza- 
tion is normalized to its equilibrium value. 



with an amplitude in the material E' Q of 2.6 x 10 8 V/m for nickel and 2.9 x 10 8 V/m for 
iron.— We take these values of the field amplitude as constant throughout the sample for the 
calculation and neglect the attenuation due to absorption in the material as well as additional 
reflection/absorption due to oxide and protection layers. We therefore overestimate the field 
amplitude present in samples used for the experimental determination of the magnetization 
dynamics. 

The optical excitation contribution alone, i. e., the second term in the dynamical equa- 
tion fll]), leads to the magnetization dynamics shown in Fig. [TJ This result should be com- 
pared to experimental values for a pulse energy density of 4mJ/cm 2 , such as those reported 



in Refs. 



10 



and !oJ where a "quenching" of the magnetization down to 40 % and 80 % of the 
equilibrium magnetization was found for nickel and iron, respectively. It is clear from Fig. [T] 
that the magnetization change computed including only the incoherent optical excitation 
at the photon energy of 1.55 eV is orders of magnitude smaller than the one observed in 
experiment. 

As a contribution to the magnetization change, the optical excitation is negligible, but it 
is still interesting to take a closer look at the carrier distributions created by the laser pulse 
because these are essentially the starting point of the momentum-resolved electron-phonon 
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Figure 2. (Color online) Spin-resolved density of states for nickel and iron (a) as well as energy- 
and spin-resolved occupation change after the optical excitation for nickel (b) and iron (c). 
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scattering dynamics. To this end, we analyze the energy- dependent occupation changes 

AN a (e,t) = J^5(e - e>i)(V^ [nj(t) - /(ej,r )] (23) 

between the dynamical distributions n-t(t) compared to the equilibrium Fermi-Dirac func- 
tions /(e^, To) that describe the carrier distributions in equilibrium at the sample tempera- 
ture T before the optical excitation. In Eq. fl23|) we separate occupation changes for minority 
and majority spins, i. e., for o = + and — , respectively, by projecting on the majority and 
minority spin contributions of the AS W wave functions using the spin-dependent weights 
(projections) 

L 

of each state ib'i where 

K 



(25) 



(^I^)l = AZ(k)Al a (k)(]\j) la + CZ(k)Al a (k){h\j)i« 
+ AZ{k)Cl c {k)Q\h)i„ + CZ(k)C^(k)(h\~h) la . 
The overlaps are given by 

(f\9)ia = / r 2 f lrT (r)~gUr)dr. (26) 
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In this paper, we always use To = 300 K as a starting point for the dynamical calculations 
to facilitate comparison with typical room-temperature measurements. Moreover, room 
temperature is still less than half the Curie temperature so that we can expect the exchange 
splitting to be not too different from its T = K value, and therefore the DFT band structure 
should be a reasonable approximation. 

Figures El^b) and E(c) show the spin- and energy-resolved occupation change, computed 
according to Eq. ( 123|) . due to optical excitation at times well after the pump pulse. It 
corresponds to the magnetization shown in Fig. [1] at 200 fs. In both materials, mainly 
minority carriers are excited. The pronounced negative and positive spikes in the minority- 
spin occupation changes are separated by the photon energy 1.55 eV and roughly coincide 
with maxima of the density of states [see Fig. E]Ja)] for the minority carriers. These maxima 
stem from the d-bands in these materials, which leads us to conclude that they play a major 
role in the optical excitation process. 

The information about the distribution after the optical excitation contained in Figs. EJ^b) 
and[2](c) allows one to draw conclusions about the maximal magnetization change achievable 
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by electron-phonon scattering in our model, as this distribution is the starting point for the 
scattering dynamics. Electron-phonon scattering is a quasi-elastic process involving a single- 
electron, i. e., there is only a small amount of energy transferred in each scattering event. Due 
to the bath assumption for the phonon system in our calculation, there are also no "secondary 
electrons" excited because such a transfer of energy to other electrons could only happen 
mediated by a phonon. We therefore expect that electron-phonon scattering will lead to a 
continuous relaxation of the excited carriers where the number of non-equilibrium electrons 
and holes decreases as they are scattered towards the Fermi energy. The demagnetization 
itself is caused by "spin- flip" scattering processes, i. e., scattering processes with different 
spin expectation values for initial and final wave functions, that occur during the relaxation 
process. Therefore the maximal demagnetization that can be caused by scattering in a fixed 
band structure occurs when all excited majority electrons and minority holes flip their spin 
while the minority electrons and majority holes do not undergo spin-flip scattering. The 
relative magnetization change in this physically rather unreasonable case is then given by 

m M = ^-2(N> + N l)tlB (2?) 

Here, \is is the Bohr magneton, and M eq and fi eq denote the equilibrium values of the 
material magnetization and of the magnetic moment per unit cell, respectively. The number 
of majority electrons iV_ e (minority holes NY) per unit cell can be obtained from integrating 
the occupation changes in Fig. [2] above (below) the Fermi energy. With that estimate we 
find a minimal relative magnetization due to electron-phonon scattering of 0.84 in nickel and 
0.94 in iron, which is a smaller demagnetization than observed in experiments. Without even 
calculating the full dynamics, we thus expect that microscopic electron-phonon scattering 
with a fixed band structure is not responsible for the pronounced drop of the magnetization 
observed in experiments. 

We next take a closer look at the photon energy dependence of the excitation process. 
Fig. [3] shows the magnetization change M^/M^ vs. the pump photon energy for a fixed 
electric field amplitude of 2 x 10 8 V/m. For most of the pump photon energies between 
0.1 and 4.5 eV the optical excitation leads to negligible demagnetization, as was discussed 
above in connection with Fig. CD for a pump photon energy of 1.55 eV. Only at a pump 
photon energy of about 0.7 eV for nickel and 2.2 eV for iron one observes a magnetization 
change of significantly more than 1 %, because these energies correspond to the exchange 
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Figure 3. (Color online) magnetization change achievable by optical excitation alone for a range 
of pump-photon energies. 

splitting of the <i-bands in these materials, so that absorption at these energies is likely to 
be associated with a change of carrier spin. The order of magnitude of our results for the 
achievable magnetization change by optical excitation alone seems to be in agreement with 
those of Zhang et al.- for the change of magnetization due to the coherent excitation by an 
optical pulse. 

B. Electron-phonon scattering 

In this section, we present results for the carrier dynamics including both optical exci- 
tation and electron-phonon scattering. We start by examining in Fig. H] the magnetization 
dynamics and the time evolution of the energy in the electronic system for the same pa- 
rameters that were used in Fig. [1] for the case of optical excitation only. Comparing the 
magnetization dynamics including electron-phonon scattering in Fig. H] to those without, 
cf. in Fig. (U one notices a demagnetization of about 3-5%. This magnetization change is 
smaller than the estimate of the previous section. Due to the scattering, the dynamics now 
also include a relaxation to equilibrium. This can be made visible by monitoring the energy 
in the electronic system [Fig. 11(b)], which nicely shows the sudden energy transfer from the 
laser pulse and a subsequent almost exponential decay with time constants of about 2 ps for 
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Figure 4. (Color online) Normalized magnetization dynamics (a) and energy difference to equilib- 
rium of the electronic system (b) after the optical excitation including electron-phonon scattering. 
Results obtained including the spin-orbit coupling contribution in the electron-phonon matrix ele- 
ment are labeled "with SOC-ME." 

nickel and 2.5 ps for iron. These time constants are significantly longer than the electron- 
phonon coupling times obtained from the analysis of experimental data for these materials 
ranging from 0.3-0.5 ps.- 1 ^ Likely, this is because we neglect other scattering mechanisms 
(such as electron-electron scattering), which open up additional scattering paths and lead to 
an overall speed-up of the relaxation process. The magnetization for nickel even rises above 
its value at equilibrium, which is understandable because there is no fundamental law that 
prohibits non-equilibrium scattering dynamics from going through intermediate states with 
an increased magnetization. Whether these are reached depends on the band structure, the 
properties of the states involved, and the initial/excitation conditions. 

When comparing the calculated magnetization dynamics to experimental results, one 
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Figure 5. (Color online) Energy- and spin-resolved occupation changes AN a at different times for 
nickel, as shown in Fig. 01 including the spin-orbit coupling in the electron-phonon matrix element. 
The representation is in analogy to the one for the optical excitation [Figs. EJJb) and[2fc)]. 

should keep in mind that our calculation neglects changes in the band structure, i.e., the 
exchange splitting, and the subsequent relaxation of these changes back into equilibrium. 
Processes associated with the a change of the exchange splitting are expected to dominate 
the dynamics after a quasi-equilibrium magnetization has been established, namely for times 
longer than about 5ps.— A meaningful comparison with experiment of the present model 
should therefore be limited to a few picoseconds, which is the dynamical time scale for which 
the different microscopic models mentioned in the introduction have been proposed. In that 
time window, we find that roughly the same results (which for clarity reasons are only shown 
for nickel in Fig. H]) are obtained if the spin-orbit term in the interaction matrix elements 
[cf. Eq. (USD] is neglected. 

To get a better understanding of the demagnetization in the present model, it is instruc- 
tive to look at the carrier distributions at different stages of the dynamics. We present only 
the results for nickel in this paper, as the carrier dynamics in iron shows similar behavior. 
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Figure 6. (Color online) Theoretical limit for the minimal magnetization achievable by a pure 
redistribution in a fixed band structure for a range of deposited energies AE. 

In Fig. [5] we plot the energy- and spin-resolved occupation change at different times. Note 
that after the end of the optical excitation at about 50 fs the excited carrier density above 
the Fermi energy (at OeV) changes only very slowly. In contrast, there is a strong change 
in the density of holes around 1.4 eV below the Fermi energy. It is the spin-flip of these 
minority holes that leads to the demagnetization of the material in the present model. The 
faster dynamics of the holes compared to the excited electrons are due to the difference in 
the density of states, cf. Fig. EJ^a), which is considerably higher below the Fermi energy than 
above, so that in this energy region there is a larger scattering phase space. In addition, 
in that energy region there are "spin hot spots," i. e., points in the Brillouin zone where 
the states are completely spin-mixed. Their presence also contributes to spin-flip scattering 
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processes.— 



C. Qualitative considerations 

As we saw from the results of the last section that electron-phonon scattering alone 
cannot explain the experimentally observed demagnetization, the next important step seems 
to be to extend the existing model to other scattering mechanisms, e.g., electron-electron 
or electron-impurity scattering. However, an argument based on energetics shows that their 
inclusion is not likely to improve the description much, if one retains the limitation that the 
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model contain only scattering, i.e., the redistribution of carriers in a fixed band structure. 
This conclusion is based on the simple observation that demagnetization in a fixed band 
structure naturally costs energy as it requires a transfer of occupation from majority states 
to minority states which are shifted up in energy by the exchange splitting. One can make 
this observation quantitative by finding the minimal magnetization that the material can 
attain given a fixed amount of deposited energy AE. This leads to a linear optimization 
problem 

min ryn?(5 2 )? (28) 

k ~ k- £ n 

with the following constraints: 

EE^ = A ^ ( 29a ) 

J2J2 n K^ E ^ + AE ( 29b ) 

k M 

Here N eq denotes the total number of carriers and E eq the total energy of the system in 
equilibrium, i.e., before the arrival of the laser pulse. As before, the contribution from 
orbital angular momentum to the total magnetization is neglected. We solve this problem 
with the ab-initio results at hand for a range of deposited energies AE, and show the results 
in Fig. |6j Note that we present the normalized magnetization, i.e., the minimum obtained 
from the solution of Eq. (|28|) divided by the equilibrium magnetization because this value 
can be readily compared to the demagnetization measured in an experiment. These values 
represent the minimal magnetization for a carrier distribution in the fixed (equilibrium) 
band structure given the deposited energy. It holds for all scattering mechanisms that could 
be creating this distribution provided that they either conserve energy (such as electron- 
electron scattering) or lead to a loss of energy by transferring it to other systems (such as 
electron-phonon scattering) . 

By comparing the experimental demagnetization with the calculated minimal magnetiza- 
tion at the amount of energy deposited in experiment one can see whether the experimental 
results can, in principle, be explained in terms of scattering alone. This comparison turns 
out to be not so easy as quite a lot of parameters (e.g. the spot size, absorption, and re- 
flectivity) are necessary for the estimate of the deposited energy from the measured laser 
intensity and some of them are known only with a considerable uncertainty. That is why 
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we chose to estimate the deposited laser energy directly from the measured magnetization 
dynamics. This is possible if one relies on two assumptions: 

1. At about 5 ps after the laser excitation the scattering processes have locally thermalized 
the material, so that the initial non-equilibrium dynamics that started has evolved 
in a quasi-equilibrium dynamics, in which the magnetization at that time can be 
characterized by the temperature dependence of the magnetization in the ferromagnet 
M(t « 5ps) = M(T) where T = T(t w 5ps). 

2. The coupling to the substrate and other losses are so weak that almost all of the energy 
deposited by the laser is still in the material at that point (t ps 5ps). However, it has 
been evenly distributed among the inner degrees of freedom. 

These assumptions are consistent with interpretations of measured data by Koopmans et 
al.— and seem to be especially well fulfilled for measurements on thin films. They can now 
be used to extract the deposited energy from the measured magnetization at 5ps, and to 
read off the corresponding achievable minimum magnetization from Fig. [7J This can be 
compared with the "quenched" magnetization reached in the same measurement. In typical 
data for nickel^ and iron^ at high intensities we find for the normalized magnetization 
after thermalization values of MNi(5ps) = 0.1 and Mp e (5ps) = 0.8. Using the equilibrium 
temperature dependence of the magnetization M(T), we conclude that the temperature after 
local thermalization is about 625 K for nickel and about 800 K for iron, respectively. As we 
assume an even distribution among the material's degrees of freedom, we can calculate the 
deposited energy as an integral over the heat capacity C P (T): 

/■T(5ps) 

AE= dTC p (T) (30) 

</300K 

which we solved using experimental data for C p (T)2i yielding Ai?(Ni) = lOOmeV/cell and 
AE(Fe) = 160 meV/cell. For these energies, Fig.[6]yields 0.26 and 0.77 as minimal achievable 
magnetizations for nickel and iron, respectively. These values should be compared to the 
experimentally observed quenched magnetizations of 0.1 for nickel and 0.7 for iron. As 
the experimentally measured minima only slightly violate the theoretical bounds, one could 
be inclined to conclude that this argument does not rule out a demagnetization on the 
basis of pure redistribution in a fixed band structure. That view changes, however, if one 
looks at the corresponding distribution functions that are necessary to attain the theoretical 
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Figure 7. (Color online) Energy- and spin-resolved occupation distributions for nickel, (a) shows 
the distributions that is necessary to attain the minimal magnetization for AE = lOOmeV/cell 
while (b) shows a typical distribution that would be created by the optical excitation (This is 
a slightly different representation of the data shown in Fig. [2(a). Here we show the occupation 
distribution which allows an easier comparison with the equilibrium distribution). 

magnetization minima. The one for nickel is shown in Fig. [7(a) and should be compared to 
the distribution that is created by pure optical excitation [Fig. UJb)]. As discussed before, 
the optically excited distribution is the starting point for all scattering processes and we 
would expect these to bring the system back to a Fermi-Dirac distribution (at a higher 
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temperature) which is also displayed in the figure. It is not at all likely that in the course 
of this process there will be an intermediate state that has a distribution that is anywhere 
close to the one shown in Fig. [TJa) for two reasons: First, for a magnetization close to the 
theoretical limit a highly "ordered" distribution is necessary, which is unlikely to be reached 
by random scattering processes. Second, the state shown in Fig. [7(a) lies very far off from the 
direct continuous transition from the distribution in Fig. [TJb) to the equilibrium distribution, 
both in terms of a simple relaxation time approximation and if one considers a quasi-elastic 
process, such as electron-phonon scattering, where we have a slow, but continuous energy 
relaxation of the exited carriers towards the Fermi energy where eventually non-equilibrium 
electrons and holes cancel out. 

Even though this argument is not a rigorous, we find it convincing enough to draw the 
conclusion that scattering dynamics in a fixed band structure cannot explain the observed 
ultrafast demagnetization. It is therefore important to include the dynamical changes in 
the band structure, i.e., the exchange splitting, in a comprehensive microscopic theory of 
ultrafast demagnetization in ferromagnets. 

IV. CONCLUSIONS 

The main objective of this paper was to analyze in detail the dynamics due to one of 
the proposed mechanisms for ultrafast demagnetization: the Elliott- Yafet process based on 
electron-phonon scattering. To this end, we carried out a numerical analysis without ad- 
justable parameters including the laser excitation and the scattering dynamics on the level 
of Boltzmann scattering integrals. We evaluated the model for the elementary ferromag- 
nets nickel and iron utilizing realistic band structures and matrix elements obtained from 
ab-initio calculations. As in previous studies,— 1 ^ we kept the band structure fixed. In this 
case, the computed demagnetization for realistic pump-laser intensities is smaller by almost 
a factor of ten than what is observed in experiments. An additional argument shows that 
this bound for the achievable magnetization "quenching" is likely to hold as well for other 
scattering mechanisms, such as electron-electron or electron-impurity scattering. We inter- 
pret our numerical results that any fully microscopic model that tries to explain ultrafast 
demagnetization by scattering dynamics really should include a dynamical change of the 
band structure, i.e., the exchange splitting. A microscopic determination of this change 
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seems more important for the understanding of the demagnetization process than studies 
focussing on Elliott- Yafet-type mechanisms. 
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Appendix: Numerical Method 

For the numerical solution of Eq. ([I]) we replace the energy delta function in the scattering 
rates, Eq. ()!]), by a gaussian of finite width to allow a Brillouin zone integration on the chosen 
grid of fc-points. The FWHM of this broadened distribution is taken to be 15 meV for the 
DFT grid at hand, and convergence of the results with respect to grid size and distribution 
width was checked. 



1. Reducing the dimensionality of the problem 

Two simplifications help to reduce the numerical effort for the solution of Eq. ([I]): 

The first is based on the fact that only the states in a limited energy range around the 

Fermi energy will experience an occupation change in the course of the dynamics. Due to 

the structure of the equilibrium distribution 

= 7? J k T — ^ 

the states far (> 200 meV for To = 300 K) above the Fermi energy (/z w Ep) will be empty 
while those far below the Fermi energy will be fully occupied. As the incoming laser pulse 
will only cause resonant transitions between occupied and unoccupied states, the occupation 
change due to the optical excitation is limited to an energy range around the Fermi energy. 
This is clearly seen in the energy resolved occupation change due to the optical excitation in 
Fig. [2j States far away from the Fermi energy (|e^ — Ep\ ^> E^) remain at their equilibrium 
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values. As the electron-phonon scattering transfers only small amounts of energy in each 
scattering process, the occupation of these states is not influenced by the following scattering 
dynamics either. That is why one can safely assume the occupation of these states to remain 
constant in time. Only states in an energy range |e- — Ep\ < SE are actually included in 
the dynamical calculation of the occupation numbers. So, for each fc-point, we only include 
the subset of bands into the dynamical calculation that fall in the chosen energy range 
SE. For the calculations including only the optical excitation in Sec. IIII Al we took SE 
to be 5eV, but for the photon energy under investigation (1.55 eV) a much smaller range 
actually suffices. We therefore reduced it to 2 eV for the calculations including scattering in 

sec. mm 

The second simplification can be made due to the crystal symmetries of the materials un- 
der investigation. These symmetries imply that the wave functions of two states at different 
fc-points which are related by a crystal symmetry operation are also connected by the same 
symmetry operation. From that one can deduce that the modulus of an electron-phonon ma- 
trix element between two states does not change if one applies a crystal symmetry operation 
on the initial and the final state. In other words, electron-phonon scattering will not break 
the symmetry of an occupation distribution that has the same symmetry as the crystal (e.g., 
the thermal distribution before optical excitation). The optical excitation could break that 
symmetry as it involves the scalar product with an external electric field. In our paper, we 
restrict ourselves to the description of the typical experimental case where the laser field is 
parallel to the material magnetization so that the symmetry of the occupation distribution 
is not broken by the optical excitation. In this case, it is not necessary to compute all occu- 
pation numbers n't. Rather all information about the occupation distribution is contained 
in any subset {kj} of fc-points which constitutes an irreducible wedge of the Brillouin zone. 



2. Numerical evaluation 

With these two simplifications we are left with a subset of occupation numbers rij which 
need to be calculated dynamically. Here j = (kj, fij) denotes a multi-index that includes the 
/c-point as well as the band index of the state. With the help of this notation, Eq. (JT]) can 
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be reformulated to yield 



at 



dnj 



at 



+ 



opt 



at 



e-p 



E(t) | £ B°fn, + £ foAJ p (1 - n.) - (1 - n.) A^n,) , 

j i 

(A.2) 

where all n^- independent quantities are contained in the constant matrices A e_p and B opt , 
which can be precomputed for the chosen set of states. Interpreting the occupation numbers 



vector n we can alternatively write this as 



d_ 

at' 



-n 



E{t) 2 B opt n + diag (ft) A c " p (l-nj- diag (l 



n) (A c - p Y n 



(A.3) 



where diag(f?) is a matrix with the occupation numbers on the diagonal. 1 denotes a vector 
with all entries set to one. In the form of Eq. (1A.3j) the differential equation is especially 
well suited for a numerical evaluation. We used a MATLAB algorithm^ 2 , for the solution. 
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